
## REPLICATION MATERIAL FOR
## HOW DOES UNCERTAINTY AFFECT VOTERS' PREFERENCES?
## BRITISH JOURNAL OF POLITICAL SCIENCE
## AUTHOR: LOVE CHRISTENSEN


#################################################
########### FILE 5: MULTILEVEL ANALYSES #########
#################################################

library(estimatr)
library(lme4)
library(optimx)
library(merTools)
library(sjstats)
library(car)
library(StrTrim)
library(DescTools)
library(ggpubr)
library(tikzDevice)


## read the cleaned data
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
load("../data/dat_uncertainty_BJPS_cleaned.RData")

## select relevant variables for reshape

dat_ml <- dat %>%
  filter(Condition != "Control") %>%
  dplyr::select(ResponseId, riskpref, partisanship, Condition, 
         mw_likert, mw_quant, mw_ideal, mw_forecast_center, mw_forecast_uncertain_magnitude,
         ct_likert, ct_quant, ct_ideal, ct_forecast_center, ct_forecast_uncertain_magnitude,
         tpp_likert, tpp_quant, tpp_ideal, tpp_forecast_center, tpp_forecast_uncertain_magnitude) %>% 
  mutate(partsender = ifelse(Condition == "Partisan", 1, 0)) %>%
  dplyr::select(-c(Condition))

names(dat_ml)[c(7, 8, 12, 13, 17, 18)] <- c("mw_fcenter", "mw_funcertain",
                            "ct_fcenter", "ct_funcertain",
                            "tpp_fcenter", "tpp_funcertain")

## respondent_policy
dat_rp <- dat_ml %>%
  gather(variable, value, -c(1:3, ncol(dat_ml))) %>%
  separate(variable, into = c("policy", "variable"), sep = "_") %>%
  spread(variable, value) 

## reverse variables so that higher outcome values should lead to higher support
## do this for MW
dat_rp <- dat_rp %>%
  mutate(fcenter = ifelse(policy == "mw", -1 * fcenter, fcenter),
         likert = ifelse(policy == "mw", -1 * likert + 7, likert),
         ideal = ifelse(policy == "mw", -1 * ideal + 7, ideal),
         quant = ifelse(policy == "mw", -1 * quant + 15, likert))

## respondent-policy-question
dat_rpq <- dat_rp %>% 
  gather(question, response, -c(1:7))

# Model 1: Policy-Respondent-Question

model_prq <- lmer(response ~ 1 + fcenter * partsender + funcertain * partsender + 
            (1 + fcenter | policy) +
            (1 + fcenter | question),
          data = dat_rpq,
          REML = FALSE, 
          lmerControl(optimizer ='optimx', optCtrl=list(method='bobyqa')))
ranef_prq <- REsim(model_prq, n.sims = 1000)[c(4:6, 10:12), -c(5)]
ranef_prq <- ranef_prq %>%
  mutate(mid = mean + 0.0341185,
         low = mid - 1.96 * sd,
         high = mid + 1.96 * sd,
         groupFctr = ifelse(groupFctr == "policy", "Policy", "Outcome"))



plot_prq_question <- ggplot(subset(ranef_prq, groupFctr == "Outcome")) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_pointrange(aes(x = factor(groupID), 
                      y = mid, 
                      ymin = low,
                      ymax = high),
                  lwd = .8, 
                  position = position_dodge(width = 1/2)) +
  coord_flip() +
  theme_minimal(base_size = 9)  +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.text.y=element_text(colour="black", size = 10),
        axis.title.y=element_text(colour="black", size = 10),
        axis.text.x=element_text(colour="black", size = 10),
        axis.title.x=element_text(colour="black", size = 10)) +
  scale_y_continuous("\n Marginal ATE Prediction Center",
                     limits = c(-0.05, 0.1), 
                     breaks = seq(-0.05, 0.1, 0.05)) +
  xlab("") +
  ggtitle("Outcome") +
  scale_x_discrete(labels = c('Idealistic','Attitude', 'Belief'))

plot_prq_policy <- ggplot(subset(ranef_prq, groupFctr == "Policy")) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_pointrange(aes(x = factor(groupID), 
                      y = mid, 
                      ymin = low,
                      ymax = high),
                  lwd = .8, 
                  position = position_dodge(width = 1/2)) +
  coord_flip() +
  theme_minimal(base_size = 9)  +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.text.y=element_text(colour="black", size = 10),
        axis.title.y=element_text(colour="black", size = 10),
        axis.text.x=element_text(colour="black", size = 10),
        axis.title.x=element_text(colour="black", size = 10)) +
  scale_y_continuous("\n Marginal ATE Prediction Center",
                     limits = c(-0.05, 0.1), 
                     breaks = seq(-0.05, 0.1, 0.05)) +
  xlab("") +
  ggtitle("Policy") +
  scale_x_discrete(labels = c('Corporate Tax','Minimum Wage', 'TPP'))


tikz(file = "../output/figures/fig_prq_results.tex", width = 9, height = 4.5)
ggarrange(plot_prq_question,
          plot_prq_policy,
          nrow = 1)
dev.off()

####################################
#####  Policy-Respondent Level #####
####################################

pr_quant <- lmer(quant ~ relevel(factor(fcenter), ref = 5) * partsender + 
            factor(funcertain) * partsender +
            (1 | policy),
          data = dat_rp)

pr_likert <- lmer(likert ~ relevel(factor(fcenter), ref = 5) * partsender + 
            factor(funcertain) * partsender +
            (1 | policy),
          data = dat_rp)


pr_ideal <- lmer(ideal ~ relevel(factor(fcenter), ref = 5) * partsender + 
            factor(funcertain) * partsender +
            (1 | policy),
          data = dat_rp)

dat_quant_exp <- data.frame(coefname = rownames(summary(pr_quant)$coef)[2:9],
                        sender = "expert",
                        variable = "PredictionCenter",
                        outcome = "Belief",
                        value = c(-4:-1, 1:4),
                        mean = summary(pr_quant)$coef[2:9, 1],
                        se = summary(pr_quant)$coef[2:9, 2]) %>%
  mutate(low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_quant_part<- data.frame(coefname = rownames(summary(pr_quant)$coef)[13:20],
                               sender = "partisan",
                               variable = "PredictionCenter",
                            outcome = "Belief",
                               value = c(-4:-1, 1:4),
                               mean = summary(pr_quant)$coef[2:9, 1] +
                              summary(pr_quant)$coef[13:20, 1])

dat_quant_part <- dat_quant_part %>%
  mutate(coefname2 = paste("re", StrTrim(coefname, pattern = ":partsender"), sep = ""))

dat_quant_part <- dat_quant_part %>%
  mutate(cov = apply(dat_quant_part, 1, 
                     function(dat) vcov(pr_quant)[dat["coefname"], dat["coefname2"]]),
         var1 = apply(dat_quant_part, 1, 
                      function(dat) vcov(pr_quant)[dat["coefname2"], dat["coefname2"]]),
         var2 = apply(dat_quant_part, 1, 
                      function(dat) vcov(pr_quant)[dat["coefname"], dat["coefname"]]),
         se = sqrt(var1 + var2 + 2 * cov),
         low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_quant <- rbind(dat_quant_exp,
                   dat_quant_part)


dat_likert_exp <- data.frame(coefname = rownames(summary(pr_likert)$coef)[2:9],
                            sender = "expert",
                            variable = "PredictionCenter",
                            outcome = "Attitude",
                            value = c(-4:-1, 1:4),
                            mean = summary(pr_likert)$coef[2:9, 1],
                            se = summary(pr_likert)$coef[2:9, 2]) %>%
  mutate(low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_likert_part<- data.frame(coefname = rownames(summary(pr_likert)$coef)[13:20],
                            sender = "partisan",
                            variable = "PredictionCenter",
                            outcome = "Attitude",
                            value = c(-4:-1, 1:4),
                            mean = summary(pr_likert)$coef[2:9, 1] +
                              summary(pr_likert)$coef[13:20, 1])

dat_likert_part <- dat_likert_part %>%
  mutate(coefname2 = paste("re", StrTrim(coefname, pattern = ":partsender"), sep = ""))

dat_likert_part <- dat_likert_part %>%
  mutate(cov = apply(dat_likert_part, 1, 
                     function(dat) vcov(pr_likert)[dat["coefname"], dat["coefname2"]]),
         var1 = apply(dat_likert_part, 1, 
                      function(dat) vcov(pr_likert)[dat["coefname2"], dat["coefname2"]]),
         var2 = apply(dat_likert_part, 1, 
                      function(dat) vcov(pr_likert)[dat["coefname"], dat["coefname"]]),
         se = sqrt(var1 + var2 + 2 * cov),
         low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_likert <- rbind(dat_likert_exp,
                   dat_likert_part)


dat_ideal_exp <- data.frame(coefname = rownames(summary(pr_ideal)$coef)[2:9],
                             sender = "expert",
                             variable = "PredictionCenter",
                             outcome = "Idealistic",
                             value = c(-4:-1, 1:4),
                             mean = summary(pr_ideal)$coef[2:9, 1],
                             se = summary(pr_ideal)$coef[2:9, 2]) %>%
  mutate(low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_ideal_part<- data.frame(coefname = rownames(summary(pr_ideal)$coef)[13:20],
                             sender = "partisan",
                             variable = "PredictionCenter",
                             outcome = "Idealistic",
                             value = c(-4:-1, 1:4),
                             mean = summary(pr_ideal)$coef[2:9, 1] +
                               summary(pr_ideal)$coef[13:20, 1])

dat_ideal_part <- dat_ideal_part %>%
  mutate(coefname2 = paste("re", StrTrim(coefname, pattern = ":partsender"), sep = ""))

dat_ideal_part <- dat_ideal_part %>%
  mutate(cov = apply(dat_ideal_part, 1, 
                     function(dat) vcov(pr_ideal)[dat["coefname"], dat["coefname2"]]),
         var1 = apply(dat_ideal_part, 1, 
                      function(dat) vcov(pr_ideal)[dat["coefname2"], dat["coefname2"]]),
         var2 = apply(dat_ideal_part, 1, 
                      function(dat) vcov(pr_ideal)[dat["coefname"], dat["coefname"]]),
         se = sqrt(var1 + var2 + 2 * cov),
         low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_ideal <- rbind(dat_ideal_exp,
                    dat_ideal_part)

dat_nonlinear_pc <- rbind(dat_ideal,
                       dat_quant,
                       dat_likert)

plot_pc <- ggplot(subset(dat_nonlinear_pc, variable == "PredictionCenter"), 
       aes(colour = outcome,
           shape = sender)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_pointrange(aes(x = factor(value), 
                      y = mean, 
                      ymin = low,
                      ymax = high),
                  lwd = .8, 
                  position = position_dodge(width = 1/2)) +
  coord_flip() +
  theme_minimal(base_size = 9)  +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.text.y=element_text(colour="black", size = 10),
        axis.title.y=element_text(colour="black", size = 10),
        axis.text.x=element_text(colour="black", size = 10),
        axis.title.x=element_text(colour="black", size = 10),
        strip.text.x = element_text(size = 10)) +
  facet_wrap(~outcome) +
 scale_y_continuous("\n Marginal ATE",
                    limits = c(-0.8, 0.8), 
                    breaks = seq(-0.8, 0.8, 0.4)) +
 xlab("") +
 ggtitle("")


dat_quant_exp <- data.frame(coefname = rownames(summary(pr_quant)$coef)[11:12],
                            sender = "expert",
                            variable = "Spread",
                            outcome = "Belief",
                            value = c(1:2),
                            mean = summary(pr_quant)$coef[11:12, 1],
                            se = summary(pr_quant)$coef[11:12, 2]) %>%
  mutate(low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_quant_part<- data.frame(coefname = rownames(summary(pr_quant)$coef)[21:22],
                            sender = "partisan",
                            variable = "Spread",
                            outcome = "Belief",
                            value = c(1:2),
                            mean = summary(pr_quant)$coef[11:12, 1] +
                              summary(pr_quant)$coef[21:22, 1])

dat_quant_part <- dat_quant_part %>%
  mutate(coefname2 = StrTrim(coefname, pattern = ":partsender"), sep = "")

dat_quant_part <- dat_quant_part %>%
  mutate(cov = apply(dat_quant_part, 1, 
                     function(dat) vcov(pr_quant)[dat["coefname"], dat["coefname2"]]),
         var1 = apply(dat_quant_part, 1, 
                      function(dat) vcov(pr_quant)[dat["coefname2"], dat["coefname2"]]),
         var2 = apply(dat_quant_part, 1, 
                      function(dat) vcov(pr_quant)[dat["coefname"], dat["coefname"]]),
         se = sqrt(var1 + var2 + 2 * cov),
         low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_quant <- rbind(dat_quant_exp,
                   dat_quant_part)

dat_ideal_exp <- data.frame(coefname = rownames(summary(pr_ideal)$coef)[11:12],
                            sender = "expert",
                            variable = "Spread",
                            outcome = "Idealistic",
                            value = c(1:2),
                            mean = summary(pr_ideal)$coef[11:12, 1],
                            se = summary(pr_ideal)$coef[11:12, 2]) %>%
  mutate(low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_ideal_part<- data.frame(coefname = rownames(summary(pr_ideal)$coef)[21:22],
                            sender = "partisan",
                            variable = "Spread",
                            outcome = "Idealistic",
                            value = c(1:2),
                            mean = summary(pr_ideal)$coef[11:12, 1] +
                              summary(pr_ideal)$coef[21:22, 1])

dat_ideal_part <- dat_ideal_part %>%
  mutate(coefname2 = StrTrim(coefname, pattern = ":partsender"), sep = "")

dat_ideal_part <- dat_ideal_part %>%
  mutate(cov = apply(dat_ideal_part, 1, 
                     function(dat) vcov(pr_ideal)[dat["coefname"], dat["coefname2"]]),
         var1 = apply(dat_ideal_part, 1, 
                      function(dat) vcov(pr_ideal)[dat["coefname2"], dat["coefname2"]]),
         var2 = apply(dat_ideal_part, 1, 
                      function(dat) vcov(pr_ideal)[dat["coefname"], dat["coefname"]]),
         se = sqrt(var1 + var2 + 2 * cov),
         low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_ideal <- rbind(dat_ideal_exp,
                   dat_ideal_part)

dat_likert_exp <- data.frame(coefname = rownames(summary(pr_likert)$coef)[11:12],
                            sender = "expert",
                            variable = "Spread",
                            outcome = "Attitude",
                            value = c(1:2),
                            mean = summary(pr_likert)$coef[11:12, 1],
                            se = summary(pr_likert)$coef[11:12, 2]) %>%
  mutate(low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_likert_part<- data.frame(coefname = rownames(summary(pr_likert)$coef)[21:22],
                            sender = "partisan",
                            variable = "Spread",
                            outcome = "Attitude",
                            value = c(1:2),
                            mean = summary(pr_likert)$coef[11:12, 1] +
                              summary(pr_likert)$coef[21:22, 1])

dat_likert_part <- dat_likert_part %>%
  mutate(coefname2 = StrTrim(coefname, pattern = ":partsender"), sep = "")

dat_likert_part <- dat_likert_part %>%
  mutate(cov = apply(dat_likert_part, 1, 
                     function(dat) vcov(pr_likert)[dat["coefname"], dat["coefname2"]]),
         var1 = apply(dat_likert_part, 1, 
                      function(dat) vcov(pr_likert)[dat["coefname2"], dat["coefname2"]]),
         var2 = apply(dat_likert_part, 1, 
                      function(dat) vcov(pr_likert)[dat["coefname"], dat["coefname"]]),
         se = sqrt(var1 + var2 + 2 * cov),
         low = mean - se * 1.96,
         high = mean + se * 1.96) %>%
  dplyr::select(variable, sender, outcome, mean, low, high, value)

dat_likert <- rbind(dat_likert_exp,
                   dat_likert_part)

dat_nonlinear_ps <- rbind(dat_ideal,
                          dat_quant,
                          dat_likert)


plot_ps <- ggplot(dat_nonlinear_ps, 
                  aes(colour = outcome,
                      shape = sender)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_pointrange(aes(x = factor(value), 
                      y = mean, 
                      ymin = low,
                      ymax = high),
                  lwd = .8, 
                  position = position_dodge(width = 1/2)) +
  coord_flip() +
  theme_minimal(base_size = 9)  +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.text.y=element_text(colour="black", size = 10),
        axis.title.y=element_text(colour="black", size = 10),
        axis.text.x=element_text(colour="black", size = 10),
        axis.title.x=element_text(colour="black", size = 10),
        strip.text.x = element_text(size = 10)) +
  facet_wrap(~outcome) +
  scale_y_continuous("\n Marginal ATE",
                     limits = c(-0.8, 0.8), 
                     breaks = seq(-0.8, 0.8, 0.4)) +
  xlab("") +
  ggtitle("")


tikz(file = "../output/figures/fig_pr_ps.tex", width = 9, height = 4.5)
plot_ps
dev.off()

tikz(file = "../output/figures/fig_pr_pc.tex", width = 9, height = 4.5)
plot_pc
dev.off()
